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Abstract 

For any real limit-n 2nth-order selfadjoint linear differential expression on [0, oo), Titchmarsh- 
Weyl matrices M(A) can be defined. Two matrices of particular interest are the matrices M_d(A) 
and Mjv(A) associated respectively with Dirichlet and Neumann boundary conditions at a:: = 0. 
These satisfy M_d(A) = — Mjv(A)~^. It is known that when these matrices have poles (which can 
only lie on the real axis) the existence of valid HELP inequalities depends on their behaviour 
in the neighbourhood of these poles. We prove a conjecture of Bennewitz and use it, together 
with a new algorithm for computing the Laurent expansion of a Titchmarsh-Weyl matrix in 
the neighbourhood of a pole, to investigate the existence of HELP inequalities for a number of 
differential equations which have so far proved awkward to analyse. 



1 Introduction 

In a recent paper |Q] a numerical algorithm was reported for the computation of the Titchmarsh-Weyl 
M matrices associated with the fourth order differential equation 

M[y] = {{py'y -{sy')y +qy = Xy. (1) 

The motivation for that work was an investigation of the HELP type integral inequahty 

^ {p\y" \^ +s\y' \^ +q\y\^)dxj < K ' ^ '^"'X \M[y]\Ux. (2) 

In this paper we turn our attention to higher order operators and inequalities. We consider operators 
of the form 

as we do not wish to be concerned with quasidifferential expressions, we assume that the coefficients 
Pj are all smooth. We assume that Pn > on (0, oo), and we assume that M is regular at x = 0. A4 
will possess selfadjoint realisations in i^[0,oo): we assume that x = oo is of limit-point (minimal 



deficiency index) type, so that these realisations depend only upon choice of boundary conditions at 
X — 0. The corresponding HELP inequality is 




< K (^^" \y\^dx^ (^j^ \M[vfdx^ . (4) 



Before we explain the connection between Titchmarsh-Weyl M-matrices and inequalities (g) and 
we give a very brief overview of M-matrices, starting with the scalar case (1x1 matrices). 
Consider a second order Sturm-Liouville equation, say 

- {py')' + qy^ Awy, (5) 

on an interval [0, oo), with a; = a regular point and x = oo a singular point of limit-point type. 
Suppose that yo denotes the solution of (^) subject to the Dirichlet conditions 

yz,(0) = 0, ra^(0) = l, 

while UN denotes the solution subject to the Neumann conditions 

2/jv(o) = -i, ra5v(o) = o. 

Then the Dirichlet m-function m£){X) is, for Q(A) 7^ 0, the unique function such that 

VNi-) + mD{X)yD{-) 

is a solution of (|^) square integrable over [0, 00) with respect to w. The Neumann m-fmiction mAr(A) 
is defined by the property that 

yoi-) ~ mN{\)yD{-) 

is a square integrable over [0, 00). The functions mu and niN are analytic functions of A on both the 
upper and lower half planes; moreover m£){X)m]\r{X) = — 1 wherever both mu and totv are defined. 

These ideas were generalized by Hinton and Shaw Q to certain Hamiltonian systems (which 
include (reformulations of) higher order Sturm-Liouville systems such as (||) - see Section 3.2 below 
for details of the fourth order case). Consider a system 

;:')-5(-;A)(:) (6) 

in which 5 is a 2n x 2n symmetric matrix given in terms of real matrices A and B by 

Six;X) = A{x) + XB{x), 

where B is positive semi-definite and, in a certain sense, positive definite on solutions of the differ- 
ential equation. The dependent variables v and u are rt-vector functions of x and A. Let 

Ud 
Vd 

be the 2n x n matrix solution of this ODE, partitioned into n x n blocks, subject to the initial 
conditions 

UDiO) = 0, VDiO) = I, 

and let 

■ Un 
Vn 
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be the 2n x n matrix solution subject to initial conditions 



Un{Q)---I: Vn{Q) = 0. 
Then the Dirichlet and Neumann M-matrices arc defined respectively by the requirements that 

Uii-):=UN{-) + MD{X)UDi-), U2{-) := Ud{-) - Mn{X)Un{-) (7) 
be square integrable with respect to B over [0, oo), for 3(A) ^ 0: 



U^{x)B{x)Ui{x)dx < +00, / U2ix)B{x)U2ix)dx < +oo. 
Jo 

The limit-point hypothesis at infinity ensures that A4n and Mo are uniquely determined 
by these conditions. 

We shall also make extensive use of the identity 

Md{X)Mn{X} = -I 

which holds wherever AId and AIn are defined. For further details see Hinton and Shaw ||^. 

Returning to the HELP type inequality we remark that this has been investigated by Russell 
[ pl , while the more general form has been studied by Dias Q . The form (||) was also investigated 
earlier in a somewhat more restricted form by Bradley and Everitt and Brodlie and Everitt B. 
In all these investigations the existence of a valid inequality, that is a finite number K in ^ or (Q), 
was shown to depend upon the behaviour of the Titchmarsh-Weyl matrix Afjv associated with (|^) or 
(||). The existence of the inequality and value of the best constant is determined by the behaviour 
of the function 

3(A^Afw(A)) (the imaginary part of A^A/at (A)) (8) 

for strictly complex values of the spectral parameter A that lie in the first and third quadrants 
of the complex plane. Indeed the existence, but not necessarily the value, of the best constant is 
determined by ^ for values such that | A |^ 0. As it is difficult to find examples of A4 matrices 
for (^ or (|^) which are known in closed form, it is of some importance to be able to investigate this 
problem numerically. 

It is further known that when lies both in the resolvent set of the realisation of A4 generated by 
Neumann boundary conditions (w(0) = in the Hamiltonian formulation, or —{py")' (0) + sy' (0) = 0, 
py"{0) = in the fourth order case) and also in the resolvent set of the realisation generated by 
Dirichlet boundary conditions (u(0) = 0, or y(0) — 0, y'{0) = in the fourth order case) then the 
inequality (Q) fails; a necessary (though not generally sufficient) condition for an inequality is that 
be a point of the spectrum of at least one of these two operators. However, when the Titchmarsh- 
Weyl matrices are meromorphic, a little more can be said on the validity of the inequality. In this 
case Dias has shown that the poles of AIn occur at the eigenvalues of the realisation of A4 
subject to Neumann conditions v{0) = and the poles of AId occur at the eigenvalues of subject 
to Dirichlet conditions u(0) ~ 0. 

If /i is a pole of an Af-matrix then fi is simple |^ and the Af-matrix has an expansion 

Af(A) = -^+ao + ... (9) 

A — /i 

where cr_i is real and is called the residue matrix of AI at fi. We shall need the concepts of Neumann 
and Dirichlet translates. In (^) we say that /i is a Neumann translate if it is an eigenvalue of 

M[y]=fiy (10) 

with Neumann condition given in the Hamiltonian form by i;(0) = 0. In a similar way we say that /i 
is a Dirichlet translate if it is an eigenvalue of (llQ), but this time with initial conditions u(0) = 0. 
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The residue matrix of A/jv associated with a Neumann translate we denote by cr at, and the residue 
matrix of AId subject to a Dirichlet translate we denote by an- 

Suppose that iin,IJ-d are Neumann and Dirichlet translates respectively. It is shown in ||^ that 
a valid inequality will be found if the differential expression is replaced by either 

Mn ^ M - ^iN 

or 

Md = M- flD 

provided either of the associated residue matrices aN, cfd is of full rank. The result for higher order 
HELP type inequalities is somewhat weaker than that for the second order problem since it is shown 
in Q for the second order classical HELP inequality that when the m function is meromorphic then 
a valid inequality exists if and only if is an eigenvalue of either the Neumann or the Dirichlet 
problem associated with that expression. In an attempt to strengthen this result for higher order 
operators Bennewitz (private communication, 1995) has proposed the conjecture that provided 

rank(crjv) + rank((T£i) ~ n (11) 

(half the order of the differential expression) then a valid inequality will be found. In this paper we 
shall prove this conjecture for the general even order HELP inequality. 

As we remarked above the existence of a valid inequality is determined by (||) as | A |^ and 
when Mat is meromorphic it must also have a pole at 0. Thus in order to investigate numerically 
the existence of a valid inequality we must compute the associated residue matrices aN^OD- In |l| 
it was noted that this was a difficult numerical problem and in section |^ we report on some new 
algorithms to solve it. 

In section ^ we apply our work to some equations to determine whether or not they are likely to 
possess associated HELP inequalities. 



2 The Bennewitz Conjecture 

In order to simplify the algebra we shall assume that at least one of the matrices Af£i(A), AfAr(A) 
has a pole at the origin A = 0. The Bennewitz Conjecture is as follows. 

Conjecture 2.1 (Bennewitz) Suppose that Mn (and hence Md) is meromorphic and that 

rank{Res{MD,0)) + rank{Res{MN ,0)) = n. 

Then there is a valid HELP inequality associated with the differential expression. 

2.1 Proof of the Bennewitz Conjecture 

In order to prove this result we shall require the following lemmas. 

Lemma 2.2 The Titchmarsh-Weyl matrices Md and Mjy are symmetric matrices and are also 
Nevanlinna functions, in the sense that the matrices 3(A'/£)(A)) and 3(A//Ar(A)) are positive definite 
for 3(A) > 0. 

For a proof of this result see Hinton and Shaw ^ . 

Lemma 2.3 A necessary and sufficient condition for the existence of a HELP inequality is that 
there exist numbers 0+ G (0,7r/2), 0_ G (tt, 37r/2), p+ > and p_ > such that 

^{-XHIn(X)) > V|A| e {0,p+), arg{X) e [e+,n/2) 

and 

'^{X'^MnW) > V|A| G (0,p_), arg{X) G [9^,3tt/2). 
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This result is proved in Dias' thesis and in a different form for the case rt = 2 in Russeh p| . 
Notes 



1. In Lemma 2^ and throughout the rest of this paper, we foUow the usual convention for matrices 
that relations of the form '> 0' and '< 0' indicate positive definiteness and negative definiteness 
respectively. 

2. We could state this lemma in an equivalent form in which there would be just one number 
p > equivalent to min(/9+,yO_). However for the proof that follows this form is marginally 
more convenient. 



Lemma 2.4 Suppose that Md or Mn has a simple pole at X — with a Laurent expansion 

1 9 

-M_i + Mo + AA/i + A^A/a + • • • . 
A 

Then all the coefficients occurring in this expansion are real symmetric matrices. 

Proof The symmetry of the coefficients follows from the symmetry of and Mo, see Hinton and 
Shaw For the rest of the proof we concentrate on Mn' the proof for Mo is similar. 
^From ^ with a; = it is clear that 

Mjv(A) = [/2(0;A), 

where ( ) is a 'square integrable' solution of the Hamiltonian system (more precisely, is a 



solution for which U2 is square integrable in the sense of Hinton and Shaw). Also, it can be shown 
from (0) that 

V2ix;X)^VDix;X)-MN{\)VN{x;X), 

whence setting a; = gives 

V2iO;X) = I. 

Thus 

MN{X) = U2Vf\0;X). (12) 



Now it is easy to see that if ( ^} ) is a square integrable solution for A = /i then ' ^^(a;, A^) 



V2ix; A^) ; " " .....6^..^.. .^.^ - ^ ...... ^ 

is a square integrable solution for X — JI. Since we are concerned with problems of limit-point type, 
the square integrable solution for any 3(A) 7^ is unique up to postmultiplication by an invertible 
constant matrix. Any such matrix cancels out upon taking the combination U2V^2~^, and hence 



U2V2-\0;t^) = U2V^-\0;-p). 



Thus from (12), 



Mat (A) = Mn{X). 

This implies that the coefficients in the Laurent expansion of Mj^ about A = are real matrices. □ 

Proof of the Bennev^ritz Conjecture Under the hypotheses of the conjecture, we can expand 
Mjj and A/jv in Laurent series about the point A — 0, 

Md{X) = X-Hl^i + Ma + XMi+0{X'^), (13) 
Mjv(A) = A-iA#_i + Afo + AA#i+0(A2), (14) 
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the expansions being valid in a neighbourhood of A = 0. By Lemma 2.4 the coefficients in these 
expansions are real symmetric matrices. We see from the hypothesis of the conjecture that 

rank(Af_i) + rank(Af_i) = n, 

so let rank(M_i) — r, rank(M_i) — n — r. We also know that M^Mn — M^Md — —I for all 
3(A) > (see Hinton and Shaw Multiplying the Laurent expansions, we obtain the following 
conditions. 

M_iAf_i = M-iM-i = 0, (15) 

M_iMo + MoM-i M-iMo + MqM-i = 0, (16) 

MqMq + Af_iAfi + MiM^i = MqMq + M^iMi + AfiAf_i = (17) 

Thus the columns of Af„i are orthogonal to the columns of A/_i. We can therefore choose an 
orthonormal basis of R", say {vi, . . . , v^}, such that 

Spaniel, . . . , Vr} = Span of the columns of Af_i, 

Span{wr+i, ■ ■ ■ ,'Vn} ~ Span of the columns of Af_i. 

Let V be the nxr matrix with columns vi, . . . ,Vr, and let V be the nx (n — r) matrix with columns 
Vr+i, ■ ■ ■ ,Vn- We make the following observations. 



• From (|15|), 

A/_it> = 0; V^M^i = 0; M^iV = 0; V'^ M_i = 0. (18) 

• Any vector v G R" can be written in the form 

v^Va + Va, (19) 

where a ^IC and a G R""*". 
Using (^ we have 

v^Mnv = a'^V'^MNVa + a^V'^MNVa 

+ 2Q^V'^MNVa, (20) 

where we have used the symmetry of AIn to simplify the last term. We now use the Laurent 
expansion ( [l^ ) together with the conditions (|l^ to simplify this expression. We observe that 

= V'^MqV + \V'^MiV + 0{\^). 

Now combining (|l|) and (^ we obtain M^iMqM^i = 0, which implies that V'^MqV = 0. Thus 

V'^MnV XV^MiV + 0{X'^). (21) 

At this stage it is not clear that the matrix V^MiV is of full rank, so the O(A^) terms might not 
be negligible. We shall show shortly that in fact V'^MiV is of full rank. The next term we need 
to examine is V'^Mj^V. This is much easier to deal with; the Laurent expansion ( p^ immediately 
gives 

V^MnV = \v'^M-iV + 0(1), (22) 
A 
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and the leading order term here is of fuU rank since the span of the columns of V is the same as the 
span of the columns of Af-i, and M_i is symmetric. Finally, we treat the term V'^MnV. Using 
and (|l|) we obtain 

V'^MnV ^V'^MqV + XV^MiV + 0{X'^). (23) 

Once more, since we know little about the ranks of the coefficients in this expansion, it is not 
clear that the O(A^) terms are negligible, a point which will have to be borne in mind later on. 
Substituting (|l]), (|2|) and (||) back into (|2^) we obtain 



v'^Mnv = (^XV'^ MiV + 0{X^)'^ a 



+ Qf^Af-if + 0(1) j a 

+ 2a^ (v'^MoV + XV^MiV + 0{X^)j a. (24) 

/,From the Nevanlinna property of Mn we know that we must have v'^Q{Mj^)v > for all 3(A) > 0, 
for all non-zero v G R". Choosing a = we see that this implies, in particular, that 

a'^^i^XV'^MiV + O{X^)ja>0 (25) 

for all non-zero a £ R'^. Now suppose that V'^MiV is not of full rank. Then we can choose a 
non-zero a such that 

a^{V^MiV)a = 0. 
Suppose now that with this choice of a we have an expansion of the form 

v'^^{Mn)v = a^(y^MpF)aS(Af ) -I- 0{^{XP+^)). 

where p > I and the coefhcient {V^ MpV)a is non-zero. Because p > 1, we know that 3(Ap) 
is not of one sign on the upper half plane. Thus v^^{Mi\[)v cannot be of one sign in the upper 
half plane, contradicting the Nevanlinna property of Afjy. We have thus established that the matrix 
V^MiV is of full rank; also, from (|25|), we have therefore established that it is positive definite: 

V^MiV > 0. (26) 

Next we choose a = 0, ci ^ in (p^): since S(i) < when 5(A) > 0, we see that V^M^iV is 
negative definite, i.e. 

V^M-iV < 0. (27) 

The two results (^) and (^7|) - together with ( p^ - imply the Bennewitz conjecture, as we shall 
show in the remainder of the proof. 

^From Lemma 2.3, we first need to show that there exists 9+ e (0,7r/2) and p+ > such that 
for |A| £ (0,p+) and arg(A) G [6'+,7r/2), the matrix 3(— A^AfAr(A)) is positive definite. 

With V = Va + Va as before, ( |2^ ) gives 

v^Q{-X'^Mn)v = a^3(-A^y'^Afiy + 0(A*))a 
+ a^3(-Ay^M_i\> + 0(A2)) a 

+ 2a^^ (^-X^V^MoV - X^V^MiV + 0(A^)) a. (28) 

Let A = pe*^ Fix a number 6*1 e (7r/3,7r/2). Then - sin(36') > for 6* £ [6*1,77/2]. Thus (|^) 
implies that there is a constant ui > and a number pi > such that for all p G (0,pi) and 
£ [9i, 7r/2] the first term in (|2l) satisfies 



7S (-A^y^MiF + (^(A*)) a > ujip^Waf. (29) 
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Similarly, fix 6*2 e (0,7r/2). Then sine* > for 6 e [6*2, tt - 02]- Thus @ implies that there is a 
constant uj2 > and a number p2 > such that for all p G (0, P2) and G [62, n — ^2] the second 
term in ( |28| ) satisfies 

d^S (-AF^Af_if + 0(A2)) d > W2/5||df . (30) 

We now deal with the last term in ( p8| ) . Clearly there exist positive constants Ci and C2 and r such 
that for aU p = |A| G (0, r) and 9 = arg(A) G (0, tt). 



2a^3 (-AV^MqV^ - X^V^MiV + 0(A*)) d < (Cip^l sin(26i)| + Cap' 



"«lll|d|| 



We bound the second part using Young's inequality: 

1 



C2p'\\a\\\\a\\<^C2p{p'\\af + p\\a\f}^ 



We also bound the first part using Young's inequality: 



Cip^l sin(20)|||a||||d|| < ici| sin(2(?)| {p^a^ + pM'} ■ 



Combining these inequalities we obtain 



|2 J ^1 3|„- foaW I ^2 4I ||.||2 J C'l I , ,| C2 2 



< liar I -fp'l sin(20)| + -fp^j + W&r I^pI sin(20)| + j ' ^^l) 

We now combine (^), ( |3C| ) and (^l]). Let 9* = max(6'i,02) £ (0, 7r/2) and let p* = min(pi, p2, r). 
Then for all p = |A| G (0, p*) and for aU 9 = arg(A) G [9*,tt/2] we have, from (||), 

„To/ \2n^. N„ ^ ^loi^COflM _ .3||^,||2 , f ^loi^COflM _ .ll/5,l|2 



z;^ 3(-A^M^)z; > - ^| sin(20)| - ^pj p^||«r + [^2 - -f\smi29)\ - ^pj pUd^ (32) 

Choosing 0+ G [6'*,7r/2) sufficiently close to 7r/2 (to make | sin(20)| small) and choosing p+ G (0,p*] 
sufficiently small, we can ensure that 

^^3(-A^M^).>^p3||af + ^p||dr (33) 

for aU p = |A| G (0,p+) and 9 = arg(A) G [6'+,7r/2]. This implies that '^{~X'^Mn{X)) is positive 
definite for all such A, which deals with the first condition in Lemma 

To verify the second condition in Lemma 2.3 we must show that there exists 0_ G (7r,37r/2) 
and p_ > such that 9(— A^MAr(A)) is negative definite for |A| G (0,p_) and arg(A) G [0_,37r/2). 
Looking back at the proof above, it is clear that there are only two changes to the reasoning. First, 
we need to replace (^9|) by a result of the form 

(^-X^V^MiV + 0{X^)j a < -wgp^Haf , (wg > constant), (34) 

to hold for aU |A| G (0,p3) and arg(A) G [tt + 6'3,37r/2] for some pa > and 6*3 G (0, 7r/2). This 
we can do because we can choose 6*3 G (0, 7r/2) such that — sin(3^^) < for all G [tt + ^3, 37r/2]. 
Secondly, we need to replace (^) by a result of the form 

d^!3(-Al>^M-iy + 0(A2)) d < -W4p||d||^ (w4 > constant), (35) 
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to hold for all |A| e (0,^04) and arg(A) G [tt + 6*4,271 - ^4] for some p4 > and 6*4 e (0,7r/2). This 
we can do because we can choose 64 G (0,7r/2) such that sin(6') < for all e [tt + 64, 2n — 94]. 
Choosing 9- > tt + max{93,94) sufficiently close to 37r/2 (to make |sin(26')| sufficiently small, as 
reasoned for (p3[)) and choosing p„ G (0, min(p3, ^4)] sufficiently small, one obtains an inequality of 
the form 

v^^{-X'M^)v < -^p^llaf - ^p\\ar (36) 
for all p = |A| e (0,/9_) and 9 = arg(A) £ [6'_,37r/2]. This implies that 5(-A2MAr(A)) is negative 



definite for all such A. Both the conditions in Lemma 2.3 have now been verified, and our proof is 
complete. □ 

2.2 A note on the converse of the Bennewitz Conjecture 

It seems appropriate to indicate here why we have been unable to prove the converse of the Bennewitz 
Conjecture: namely, that if 

rank(Res(Af£), 0)) + rank(Res(Afjv, 0)) = n - d < n 
then there is no HELP inequality associated with the differential operator. The key to proving 



such a converse would be Lemma 2.3, which is an if-and-only-if result. Following the notation and 
proof of the previous section, suppose that V is an n x r matrix whose columns are r orthonormal 
vectors spanning the column space of the matrix M^i and let y be an 71 x (n — r) matrix whose 
first n — r — d columns are orthonormal vectors spanning the column space of and whose last 
d column vectors are chosen so that the columns of V are orthonormal. Then one can show that 
the columns of V and of V form an orthonormal basis of R", as before; equations ( [l9| ) and ( ^ ) can 
be shown still to hold. Looking for a failure in the first condition of Lemma ^.3| , we seek a vector v 
such that for all 9+ G (0,7r/2) and p+ > 0, the inequality 

v'^'i{-X'^MN{X))v > (37) 

fails for some arg(A) G [9^, 7r/2) and |A| G (0, p+). It seems reasonable to look for such a vector v in 

that part of R" which is not spanned by the columns of M_i and Af_i: to this end we must have 

a = in (p^, and we must also choose a such that M^iVa — 0. With these two conditions, (^J) 
becomes 

v'^'^{Mn)v = a^'^iV^MoV + XV^MiV + . . .)a, (38) 

while similarly 

v^Q{-X'^Mn)v = a^'^i-X'^V^MoV - X^V^MiV + . . .)a. (39) 

Recall that all quantities in these equations, other than A, are real. Let A — re'^. Then 3'(— A^) = 
— sin(26'), which is negative for 9 = 7r/2 — e, for small e. Thus combining ( ^7| ) and (^), a necessary 
condition for a HELP inequality to hold is that 

a^V^MoVa = 0, (40) 

whence (|8|) and ( p9| ) become 

w^3(MAr)w = a'^3(AF'^MiF + ...)«, (41) 

^'^^(-A^AfAr)!; = a^'^i-X^V'^MiV + . . .)a. (42) 

The Nevanlinna condition 3(AfAr) > for 3(A) > implies that the leading term on the right hand 
side of ( pl[ ) is strictly positive; the leading term on the right hand side of (^) is then strictly positive 
for arg(A) G (7r/3, 7r/2], which certainly does not preclude the existence of a HELP inequality: indeed, 
if it were true for all v and not just those outside the column span of A/_i and A/_i, it would say 
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that a HELP inequality definitely held. This suggests that we ought to try to prove that (|4^) must 
fail, but we have so far been unable to do this. 

(Note finally that no new information is obtained by looking for a failure in the second condition 
of Lemma p.3[ ). 



3 Computing Mn{X) and Md{\) near a pole by a change of 
variables 

We now turn our attention to the problem of computing the residue matrices of Mjv and Mjj near 
a pole. We shall assume once more that the pole is at A = 0. Also, since the numerics are the same 
for Md as they are for Mat, we shall consider a more unified problem: that of computing the residue 
matrix of an arbitrary Titchmarsh-Weyl matrix M(A) having a pole at A = 0. 

We describe the solution of the problem in four steps. In the first of these, we define a new 
matrix 4* and show that it is well-behaved near the pole of M. In the second, we explain how \1/ 
can be computed by integrating an initial value problem. In the third part we explain a simple 
extrapolation procedure which we use to determine the Taylor expansion of \1/ near the pole of M; 
finally we show how the Laurent expansion of Af may be recovered from the Taylor expansion of 



3.1 The transformation to 'safe' variables — the matrix \1/ 

Suppose that the Titchmarsh-Weyl matrix Af(A) has a pole at A = 0. We know that such a pole 
must be simple but we also know that any attempt to compute M directly, by the methods 
we described in our previous work Q, is likely to yield inaccurate results when |A| is small. In an 
attempt to circumvent this difficulty we shall define a new variable \l/ by 

* (a/ + M-i)-\ (43) 

where a is a complex constant to be chosen for convenience. 

The reason for removing the singularity in this way, rather than by using \Ad{X) as a new 
variable, lies in the need to approximate whichever variable is chosen by solving an initial value 
problem. For ^ the resulting Riccati-type ODE ( psf ) is not singular as A ^ 0. For AM(A), on the 
other hand, the corresponding Riccati equation has a A^^ singularity occurring in the quadratic 
term on the right hand side. 

In order to show that a may be chosen so that 5* has a removable singularity at A = we need 
to consider two different cases separately. The first is the case where is well-behaved at A = 0; 
the second is the case where A/^^ also has a pole at A = 0. The second of these two cases seemed, 
initially, the more pathological, since it includes the case in which the Sturm-Liouville problem with 
Neumann boundary conditions shares an eigenvalue with the same problem with Dirichlet boundary 



conditions: however Lemma 4.1 below gives a whole class of problems for which this always happens. 



Lemma 3.1 Suppose that M ^ has a removable singularity at X = Q. Then there exists a choice of 
a such that ^ has a removable singularity at X — 0. 

Proof Let M denote liniA^o ^^^^(•^); this exists by hypothesis. Clearly 5* will have a removable 
singularity at A = if and only if —a is not an eigenvalue of M. Thus the result is proved. □ 

Lemma 3.2 Suppose that Af^^(A) and A/(A) both have poles at X = 0. Then ^I^ has a removable 
singularity at X = for any a with 3(a) ^ 0, and for all but finitely many real a. Also, in the case 
n = 2, the matrix ^' has a removable singularity at X — for any non-zero a. 
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bi,i 61,2 

62,1 &2,2 



Proof We examine the 2x2 case first. For any 2x2 matrix 

B = 

let B^ denote the matrix of minors of B, so that 

B^. 



&2,2 — ^1,2 
-^2,1 61,1 

and BB^ B'^B = detS. With this notation it is clear that 

a/+(Af-i)^ 
~ det(a/ + M-i)' 

Expanding the determinant, we get 



+ atrace(M-i) + det(M-i) 



We also know that 



detM 

Since M has a simple pole, so does (n.b. this step does not generalise to the case of 
n X n matrices). Also, M^^ has a simple pole by hypothesis. Returning to (|44|), we can see that 
al + {M~^)^ has a simple pole. Thus to get a removable singularity for for any non-zero a, it 
suffices to show that trace(Af~-'^) has a pole. Let 

M-\X) = + Ala + AA/i + 0{X'^); 

A 

the only way that trace(M~^) can fail to have a pole is if trace(M_i) — 0. If this happens then 
Af_i must be of the form 

Pi' 



7 -P 

recalling that M_i is a real non-zero matrix this implies that Af_i is of full rank, and hence M must 
have a zero rather than a pole when A = 0. This contradiction proves that trace(M^^) has a pole, 
and hence that ^ has a removable singularity. 

We now turn to the case n > 2. We know that M and M^^ are analytic functions of A with 
singularities at A = 0. Let Ati(A), . . . , /^n(A) be the eigenvalues of AM~^(A). As AM~^(A) is analytic, 
and symmetric in the sense of Kato p. 120], we know from the remark at the bottom of page 



121 in 1 10 that M ^ has an analytic Schur decomposition of the form 



M-\X)R{X) = R{X)\d{X) 
A 

on a punctured neighbourhood of A — 0. Here 

D = diag(^i, . . . ,/i„), 

where /xi,...,/x„ are analytic at A = 0, while the matrix R{X) is analytic at A = and is real 
orthogonal {R^^ — R^) for all sufficiently small real A. This orthogonality of R for real A means 
that R^^ is also analytic at A = 0. To see this, observe that the only type of singularity which R~^ 
could have would be a pole. A pole would cause i?~^(A) to blow up as A approached zero through 
real values, contradicting the regularity of R by the orthogonality R^^ = R^ for real A. 
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The Schur decomposition of is clearly 

^'(A) = R{X){X-^D{X) + aI)-^R-\X), (45) 
and the eigenvalues of ^' are clearly 

(a + ■■■,{a + jf^ny^- 

These are all analytic functions of A. If j is such that ^ at A = then the corresponding 
eigenvalue of 5* clearly has a zero at A = 0. If i is such that fij has a zero of order at least 2 at 
A = 0, then the corresponding eigenvalue of ^ has a removable singularity at A = provided a ^ 0. 
If j is such that fXj has a simple zero at A = then the corresponding eigenvalue of 5* will have a 
removable singularity at A = for all but one value of a. In particular, since fij is real-valued for 
real A, the corresponding eigenvalue of has a removable singularity for 3(a) ^ 0. Whenever the 
eigenvalues of ^' are analytic, so is ^' itself, from the Schur decomposition (Esh in which R{X) and 
are analytic at A = 0. This completes the proof. 

Remark The proof for n > 2 can be extended to show that under the hypothesis of the Bcnncwitz 
conjecture ~ namely, that the ranks of the residue matrices of M and sum to n - the matrix 

has a removable singularity for any non-zero a. In other words, the result of the case n = 2 is 
recovered in this special case. □ 



3.2 The initial value problem for \E' 

In the rest of this section we shall consider the case n — 2: that of the fourth order Sturm-Liouville 
problem. We start by recalling the method proposed for the computation of the matrix AI in 
An interval [0, X] is chosen, with X suitably large; the fourth order Sturm-Liouville equation is cast 
in the form 



Jz' = 5*2:, 



where J is the symplectic matrix 



J = 



S is the symmetric matrix 



S 



V 

and z is the vector of quasi-derivatives 



/ 





-1 

















-1 


1 













V 


1 







/ 


Xw — 


q 
















— s 


1 










1 



















1/p 



(46) 



V 



y 
y' 

-{py"y 
py" 



sy 



J 



Then we consider the matrix initial value problem consisting of the differential equation 

JZ' = SZ 



12 



and the initial condition 



Z{X) = 



r 


M 








1 





V 


1 / 



Let Z{x) denote the solution of this problem, where Q < x < X . We partition Z as 

-<^' - ( ) ■ 

and form the corresponding initial value problem for the matrix UV^^. We solve this initial value 
problem, starting from x — X, to find UV^^{0). Our approximation to M is then given by 

M {UV-'^{0))-\ (47) 

(The formula would be exact if we had X — +oo.) If we replace M in (^ ) by the expression on the 
right hand side of ( ^7| ) then we get 

5- « {al + UV-\0))-\ 

Clearly, then, the process of approximating can be reduced to that of deriving an initial value 
problem for the matrix 

r(a;) := (al + UV'\x))-\ 

The initial condition is obvious: r(X) — a~^I. The differential equation is also quite straightforward 
to derive. If the matrix S in (^) is partitioned as 

S ~ ( ^^'^ ^^'^ I 

then it turns out that 

r' = - {TS2.i{i - ar) + (/ - ar)5i,2r + (/ - ar)Sis{i - ar) + r52,2r} . (48) 

We solve this equation using the NAG routine D02QGF, which allows reverse communication for 
evaluation of the right hand side of the differential equation: this helps to keep the programme 
structure simple when the right hand side is complicated. D02QGF is a variable-order, variable-step 
Adams code and is therefore able to cope with mild stiffness. In practice, we noted in Q that 
stiffness is not usually a problem unless X has been chosen much larger than necessary. 

3.3 Determining the Taylor expansion of ^ 

Determining an approximate Taylor expansion of an analytic function from numerical values of the 
function is not easy. The number of coefhcients in the expansion which can be computed reliably 
depends on the accuracy with which the function values can be computed, on the rate of decay of 
the Taylor coefficients as one proceeds up the series and, ultimately, on the precision of the machine 
arithmetic. 

Our problem is slightly compounded by the fact that we have a function ^E* with a removable 
singularity at the point around which we wish to expand it (A — 0). We cannot compute ^'(0); 
indeed we cannot compute ^'(A) for any A with zero imaginary part. Our approach has been to 
compute \l/ at a sequence of points 

X=-^, j=0,l,2,..., 
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where /x is a fixed complex number with ^ 0, and solve a Vandermonde system (using the 

algorithm of Bjorck and Pereyra Q) to obtain approximations to the Taylor coefficients. 

We shall now consider how the different sources of error and the ill-conditioning of the Van- 
dermonde system will affect the approximations to the Taylor coefficients which we obtain. For 
simplicity we shall set aside our matrix-valued function 4* and consider a complex-valued function 
/ given by a Taylor expansion 



f{z) = ao + aiz + 022^ + ■ • • + a,nz"^ + a,n+iz"^~^^ + 
Clearly the following system of equations holds: 



(49) 



ao 
ao 



ao + aifj, + a2/i 

2 



• amfJ- 



a2 ' 



m+l 

-iZ — arn+2Z 



ao + ai§r + a2 (^) 
where the function g is given by 

.9(2) = f{z) - a„ 

A Vandermonde matrix is an (rn -I- 1) x (m -I- 1) matrix of the form 



V = 



= .9(m) 
= .9W2) 

- 5(m/2") 



(50) 



/ 1 


1 


1 • 


1 


\ 




ai 








2 












al 


0^1 ■ 


.. Q,2 





\ 



where ao, ■ ■ ■ , a,yi are distinct complex numbers. If we let a = (aq, aifi, 02/^^, . . . , amlJ'™)'^ and 
g = {g{^/2"'^), g{^/2™~^), . . . ,g{ii))'^ then we can cast our system as a dual Vandermonde problem 
of the form 

F^a - g, 

where ao = 1/2"*, ai = l/2™~^,. . . ,am = 1- Since the aj are positive real numbers arranged in 
ascending order, the error analysis of Higham |^ is now applicable to the solution of this system 
by the Bjorck-Pereyra algorithm. In particular this error analysis shows, when we know g 'exactly', 
that the Bjorck-Pereyra algorithm introduces essentially no more error into a than is already implied 
by the storage of g in machine arithmetic. Since we do not have any of the special sign properties 
on the vector g which would make for a better error estimate, this suggests that the contribution of 
the machine precision e to the error (measured in the norm || • ||i) will be a term of the order 



(51) 



where C is independent of m ||^. With a machine precision of 10^^^ this suggests that, on grounds 
of roundoff alone, it will not be possible to obtain reasonable accuracy in the vector a for m much 
greater than 6. This was borne out in the experiments which we conducted. Of course we could use 
II • lloo to measure the error, instead of || • ||i; however this would make no difference, as the ratio 
||x||i/||x||oo is never greater than m for any non-zero m- vector x, and we have already seen that m 
must be quite small. 

We now turn to the contribution to the error arising from g: we do not know g exactly because 
we do not know g exactly. We must approximate g hy f. This entails an error 



9{z) - .f{z) = -am+i^" 



-2Z 



(52) 
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plus a further term arising from the numerical integration of the differential equation satisfied by /. 
process where the accuracy not is much that the error alternation Neglecting the integration error 
for the moment, we observe that ao,. . . ,a„ from (|50|), the coefficients ao,. • • jCim are the interpolation 
coefficients for the function g at the points fi, /i/2,. . . ,/i2~™. As interpolation is a linear process, 

Intcrpolant of <? = Interpolant of / + Interpolant of {g — /) 

oo 

= Interpolant of / — am^pjlnterpolant of 2'"+^}, 
p=i 

the second equality coming from (p2|). Now let 



(53) 



rn / \ ^ 



denote the mth degree polynomial intcrpolant to the function z ^ z'^+p at the points 1,2 ^ , . . . , 
2~™. By (53) the change in the computed value of a^/i'' caused by approximating g by / is given 

by 

00 



p=l 



(p) 

The coefficients cl can be computed explicitly using the Bjorck-Pereyra algorithm. The solution 
is given by 



ci"^ = E (-l)""''^^'/3m-.~,(l, 2-\ . . . , 2-'"+i)Ai"+P f ^ 

j=0 ^ 

where Pk{yo, ■ • • 7 Vm-i) denotes the sum of all products of k distinct elements of the set {yo, ■ ■ • , Hm-i}- 
We use this formula to get a bound on clJ'\ First, we consider /3fe(l, 2~^, . . . , 2"™+^). This is a sum 



m+p-j T-ij-1 



(54) 



of 



m 
k 



terms, of which the greatest is 1.2^1 . . . 2-''+^ = 2-'=('=-i)/2 Thus 



f^m—u—j (1; 2 , 

Next we tackle the term 



-m+l 



)< 



J-1 

Y[ii- 2™+''-''; 

r=0 



m 

m — V — j 



(55) 



Similarly, 



11(1-2'') 



'r=0 
r=0 

^ 2(™+p-0-i)/2)i 



r=l 

= 2^0+i)/2 J]^(l_2-'-) 

> 2^0+i)/2 -Q exp(-2.2-'^) 

r=l 

> 2^0+i)/2exp(-2). 



(56) 



(57) 



15 



Substituting (|5|), (|6|) and (|^) back into ^) we obtain the estimate 



(58) 

which simpHfies to give 

|clf)| < E ( - ■ ) |Mr+Pexp(2)2-(™-'')'/2+(™-'')(2j+i)/2-J-0-+i)/2. (59) 

By the change of dummy summation variable j ^ m ^ v — k, this gives 

\c^J'^\ ^ E ( T ) lMr+^exp(2)2-'=(^-i)/^. (60) 

For our purposes, it suffices to make a very blunt estimate at this stage: throw away the powers of 
2, and extend the summation up to fc = m. Since 



E 



k=0 

we get the bound 



771 

k 



|c(f)|<|Ain2/ir"exp(2)<exp(2)|2Ai| 



Thus, neglecting the integration error, the change in the computed vahie of ai^fx^ caused by replacing 
g by f will be bounded by 

oo 

exp(2)El«m+pl|2Mr+^- (61) 
p=i 

Provided 2/i is strictly within the radius of convergence of the power series, this will tend to zero as 
m tends to infinity; moreover, if Um+i ^ 0, the leading order term for small fi will be 

\a^+i\exp{2)\2fir+\ 

In order to assess the likely order of magnitude for a suitable value of fi it seems reasonable to ask 
that this error term be of the same order of magnitude as the error arising from roundoff. The 
constant C in ( ^l| ) is independent of m so we can neglect it; the coefficient am+i we obviously do 
not know, but if we assume that it is 0(1) then we obtain 

m.2"'e « \2nr+\ 

We have already seen that with e = 10~^^, the choice of to = 6 is likely to be the biggest possible. 
In some sense, this modest value of to justifies the assumption that Qm+i is an 0(1) quantity. It 
also gives to. 2™ e « 4 x 10^^, which suggests /i « 0.1. With to = 5, on the other hand, we have 
TO. 2™% « 1.7 X 10"^, which suggests fj, « 0.025. 

Based on these observations we devised the following algorithm for computing the first k + 1 
Taylor coefficients of the matrix 5*, where k < 5. 

1. Make the tolerance TOL for the computation of as small as possible within the constraints of 
reasonable run-times. This depends on the machine at one's disposal. 



2. Start with « 0.025 (say) and n = k + 1. 
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3. Compute approximations to ao, . . . , am', from these extract approximations to ao, . . . , Uk, the 

coefficients of interest. 

4. If m < 7, increase the value of m and compute new coefficient approximations. 

5. While the approximations seem to be improving and m < 7, keep increasing m and computing 

new approximations. 

6. When the sequence of coefficient approximations appears to start to diverge stop increasing m. 

Discard the latest (starting-to-diverge) approximations. 

7. Now regard m as fixed and start to double fi. Follow the same process as above, doubling fi while 

this seems to improve the values of ao, . . . , ak- Stop either when the user's target accuracy 
is achieved, or when the approximations seem to start to diverge, or when a doubling of fi 
would give > 0.5. Return a warning flag {IFAIL ~ 2) if the target accuracy has not been 
reached. 

The error due to integration is never explicitly controlled in this process, though step 6 should 
ensure that m is never taken large enough to magnify the integration error to an unacceptable level. 

Typically, for computing the Taylor coefficients of our matrix 5*, we might start with /i = 0.025i, 
and compute using an initial value solver with TOL ~ 10^^^. If 5* has a Taylor expansion 

^'(A) = *o + A^"! + A^*2 + A^*3 + A'^*4 H , 

we usually find that only the coefficients 5*0, . . • , ^3 can be computed with an accuracy of 10^** or 
better; the accuracy of ^3 might be 10^'*, of ^'2 about 10~^, of about 10~*, while ^'g might 
have an accuracy of 10~^°, achieved with m — 6 or m — 7 and with a value of /j, of about 0.2i. 
The precise details depend on the problem in question. The deterioration in the accuracy of the 
coefficients as one proceeds up the series need not be a problem if one intends to use them simply 
to compute values of ^'(A) by a truncated Taylor expansion for small values of |A|. 

3.4 Recovering the residue of M from \1/ 

We suppose that the first few terms of the Taylor expansion of M have been determined: 

*(A) = *o + A^-i + A2*2 + A^*3 + • • • (62) 
We want to determine the first few coefficients in the Laurent expansion of M: 

M(A) = \m_i + Mo + AAfi + \^M2 + ■■■ (63) 
A 

As we shall see, from the first n terms in the Taylor expansion of 5* we can determine at most the 
first n — 1 terms in the Laurent expansion of M. Thus, although we have eliminated the problems 
associated with trying to compute M near the pole, we have paid a price in terms of having to 
compute more Taylor coefficients than we get repaid in Laurent coefficients. 
Equation ( ^3|) may be rearranged to yield 

Af-i = - al 

d^tM""^ 
*^-adet(^')/ 
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Inverting both sides, 

M = det(1')(«''^ -adet(1')/)~i 

, ^-adetf*)/ 

= detf^*^ — 

^ ^det(*-adet(«')/) 

Now we expand the determinant in the denominator to get 

det(^' - adet(^')/) = det(*) {l - atrace(*) + det(*)} . 
Thus we obtain the formula 

^f-adetf*)/ , , 

1 - Q;trace(*) + det(^') ■ 

We shall obtain the Laurent series for M by Taylor expansion of the numerator and denominator in 
(p4[). From (^2|) we obtain the expansions 

trace(5') = trace(5'o) + Atrace('I'i) + A^trace(*2) + A^trace(v[>3) H , (65) 

det(5') = det *o + Atrace(^'o *i) 

+ A2(det(*i) +trace(*^4'2)) 

+ A^(trace(*f ^'2) + trace(^'o ^-s)) 

+ A'*(trace(*^*4) + trace(«'f ^3) + det(*2)) 

+ ■ • • 

We know that the denominator in (joj) must have a zero at A = because M has a pole at A = 
by hypothesis: thus 

1 - Q;trace(vl/) + det(*) = A(a^trace(vl/^vl/^) - atrace(v['i)) 

+ A^(q;^ det(*i) + Q;^trace(*o *2) - Q;trace(*2)) 

+ A^(a^trace(^'o *3) + a^trace(*f ^'2) - Q;trace(*3)) 

+ A''(a^trace(*o *4) + a^trace(*f ^fa) + det(*2) - atrace(*4)) 

+ ••• (66) 

For brevity we shall write this expression in the form 

1 - atracc(*) + dct(^') = Aai + A^a2 + A^aa + A'*a4 H , (67) 

where the coefRcients ai , 02 , 03 and 04 are evident by comparing ( p6| ) and ( |67| ) . We shall also write 

^ -a det(«')/ ^ Ao + XAi + A^ A2 + A''' A3 + A** A4 + • • ■ , (68) 

where 





= ^0 


- adet(^'o)/. 


Ai 


= *i 


- Q;trace(5';f^'i)/, 


A2 


= *2 


- a (det(^'i) + trace(^'Q *2)) I, 


A3 


= *3 


- a (trace(*jf ^-3) + trace(^'f *2)) I, 


Ai 


= ^4 


- a (trace(*^^'4) + trace(^'f ^3) + det(*2)) I 
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If ai 7^ 0, then we may combine (|64[), (67) and ( |68| ) to gives us the Laurent expansion 
M(A) 



+ A 



ai 



ai 



- - ^^1 + ( - - - 



(2aia3 — 02)02 04 
1. 



^ Mo 



(69) 



Since M has at worst a simple pole, when oi = then = 0. In this case the Laurent expansion 
of M becomes 



M(A) 



+ A 



-1^ 

0.2 



0.3 



02 
A2 



04 



(70) 



This happens when ^'(0) = which happens when the residue matrix Res(Af, 0) has full rank, 
giving a zero at A = 0. The code checks that the value of oi, and uses (|69| ) if |oi| > TOL, 

(|o|) if |oi| < TOL, where TOL is the tolerance used by D02QGF in the computation of ^. The 
case oi = 02 = cannot arise with 2x2 matrices: for 02 = necessarily implies Ai — 0. Since 
^f(O) = i/, we see that det *(0) ^ 0; this imphes that Af^^ = (^-^ - det(*)/)/ det(^') has a double 
zero at A = 0, implying that M has a double pole (or worse). This is not possible, as both M and 
have, at worst, simple poles. 



4 Numerical Experiments 

Our primary objective in these experiments was to compute the residues of Titchmarsh-Weyl ma- 
trices for a number of fourth order Sturm-Liouville equations and to use these, together with the 
Bennewitz Conjecture, to decide whether HELP inequalities hold for these equations. 

Before listing our example problems, we mention the following useful result. To avoid complicated 
conditions on quasiderivatives we state the result for smooth coefficients in the differential operator. 



Lemma 4.1 Suppose that C is a fourth order differential operator on the domain of functions f in 
L^[0, 00) which are four times continuously differentiable and are such that Lf G L^[0, 00). Suppose 
that C has the form L — P' where i is a second order operator 

i{y){x) = ~y"ix)+qix)yix) 

and q is twice continuously differentiable with q{0) ~ 0, q'{0) — 0. Suppose also that C is strong 
limit-point at infinity. Then any eigenfunction y of C subject to Neumann boundary conditions 
y"'{0) = = 2/"(0) which is not in the null-space of £ generates an eigenfunction z = £y of C subject 
to Dirichlet conditions z(0) = = z'{0). 

Proof Suppose y is as described, so that Cy = Xy for some real A, and let z = £y. Because y is 
not in the null-space of £, z is non-trivial. Also z{Q) = —y"{0) q{0)y{0) — because y"{0) = 
and q{0) = 0, and z'{0) = -y"'(0) -I- g(0)?;'(0) + g'(0)y(0) = since y"'(0) = and q{0) = q'{0) 
— 0. Thus z satisfies the Dirichlet boundary conditions. Clearly Cz = I'^ly — £Cy — £\y — X£y 
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= Xz, so z satisfies tlie differential equation Cz = \z: as the coefficient q is twice continuously 
differentiable, this makes it easy to see that z is four times continuously differentiable. Finally, z is 
square integrable. This follows because 

{z,z) = {iy.iy) = {ey,y) ^ \{v,y), 

the penultimate equality using the fact that l"^ — C is strong limit-point at infinity together with 
the fact that ly{Q) = (^?;)'(0) =0. □ 

Equation 1 The differential equation 

y''") - {s{x)y')' + q{x)y ^ Xy 
on the interval [0, c»), with coefficients 

Sx^ix^ -3x^ -5) 



s{x) 

q{x) = 



{X^ + 1)2 

4[4a;i2 _ 24xi° - 7x^ + 96x^ + AQx^ - eOa;^ - 15] 



(a;2 + 1)4 

This equation is strong limit-point at infinity. Imposing a Dirichlet boundary condition y(0) = 
?/'(0) = reveals that this problem was carefully crafted so that A = would be an eigenvalue 
of multiplicity 2; the reader may check that 

y{x) = {x^ + a;'*)e"^^ y{x) = {x^ + a;^)e"^' 

are the corresponding eigenfunctions. We arranged this because we suspected that it would 
result in a problem for which we would have 

rank(Mi)) = 2, rank(AfAr) 0, (71) 

although in fact we do not know of any result which would guarantee this. The numerical 
results in Table |l| suggest that (|7l] ) is indeed true (the determinant of Mn vanishes nowhere) . 
If it is, then the hypotheses of the Bennewitz conjecture are satisfied and we have an equation 
for which a HELP inequality holds. We should mention that there seems to be a dearth of 
fourth order examples with multiple eigenvalues in the literature: indeed, we could not find 
any 

Notice that we carried out the computations for two different values of a, to provide an 
additional check on our results; we also quote what the code thinks is the imaginary part 
of the residue matrix. This ought to be zero, so it provides an indication of the error. We 
also quote the error indicator returned by the code: this is reassuringly of the same order of 
magnitude as the imaginary part of the computed residue matrix. 

Equation 2 The differential equation 

y(™) - 2{xy)' + [x^ - 2)y = Xy 
on the interval [0, oo). This equation is strong limit-point at infinity; its differential operator is 



the formal square of the second order operator ty = —y" + x'^y. By Lemma 4.1 we know that 
all but at most one of the Neumann eigenvalues will be Dirichlet eigenvalues: in fact if Dirichlet 
boundary conditions j/(0) = = y'(0) are imposed, then the eigenvalues are Xk = 16(fc -I- 1)^, 
while if Neumann boundary conditions y"'{0) = = y"(0) are imposed, then the eigenvalues 
are Xk = I6k^. Thus at each of the points A — 16(fc + 1)^, fc = 0, 1, . . ., both the Titchmarsh- 
Weyl matrices Mjy and AIo will have poles. Given that Aljy = — A/t^^ it is clear that this 
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Using a — 1, truncating [0, cx)) to [0, 100]: 


Res(Mr,, A = 


. / -384.83 -167.96 \ . / 4.2 x 10"^ -1.7 x 10"^ \ 
^"1^-167.96 -79.29 j ^ * 1^ -1.7 x lO'^ 6.6 x lO^^ J 
Code error estimate (sup norm): 3.4 x 10^^ 


Determinant of residue matrix: 2303.17 + iO.0008 




Value of |ai|: 6.3 x 10"^". 




Integration tolerance for ^P: 10"^ 


Using a = 1 + i, truncating [0, oo) to [0, 10]: 


Res(Af_D, A 


. / -384.79 -167.94 \ ■ ( 2.5 x 10"^ 1.0 x 10"^ \ 
~ ~ -167.94 -79.28 y 1.0 x 10"^ 4.1 x 10"'' j 
Code error estimate (sup norm): 2.1 x 10^'^ 


Determinant of residue matrix: 2303.18 — «0.02 




Value of |ai|: -1.7 x lO^^^. 




Integration tolerance for 'J: 10^^ 



Table 1: Results for Equation 1 



Using a= 1, truncating [0, oo) to [0,20]: 


Res(M£i, A 


. / -82.62549 -40.74366 \ ( -7.3 x 10-« -3.6 x 10" 
^ ^^1^-40.74366 -20.09121 y -3.6 X 10^6 -1.8x10" 
Code error estimate (sup norm): 3.5 x 10~^ 




Determinant of residue matrix: —4.5 x 10 ^■^ — i2.3 x 10 




Value of |ai|: 9.7 x 10"^. 






Integration tolerance for ^: 10~^ 




Res(M7v, A 


. / -82.62548 -40.74366 \ ./ -1.8 x 10"'^ -8.7 x 10" 
~ ^ ~ -40.74366 -20.09121 ) \ -8.7 x lO"*^ -4.3 x 10" 
Code error estimate (sup norm): 4.7 x 10"^ 




Determinant of residue matrix: 2.3 x 10 ^■^ — i5.7 x 10 




Value of |ai|: -9.7 x 10"^. 






Integration tolerance for ^f: 10~^ 




Table 2: Results for Equation 2 



means that at each of these poles, the residue of Mn and the residue of Md will both be of 
rank 1. Together with the now-proved Bennewitz Conjecture this means that there are HELP 
inequalities (|^) associated with this differential equation, provided the operator A4 is defined 
by 

My = y^™) - 2{x^y'y + {x^ - 2)y - 16(fc + i fy, 

where k is some non-negative integer. This is a result which Diaz conjectured in his thesis |^ 
but was unable to prove. 

Consulting the numerical results in Table § we see that our code obtains approximations to 
the residue matrices which are as near to rank 1 as one could expect: they are matrices whose 
elements are not small but whose determinants are O(10~^^). 

Equation 3 The differential equation 

_ 2(e-y')' + (e'" - e^)y = Ay (72) 

on the interval [0, oo). The differential operator here is the formal square of the second order 
operator £y = —y" + exp(a;)y. With Dirichlet conditions y(0) = = y'(0) the eigenvalues 
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Using a = 1, truncating [0, oo 


) to [0,10]: 




TJ i^f \ ^^^^^^^^^ / "297.2883 -110.8968 \ 
Res(M,,, A = 35.560604) = _n0.8968 -41.3676 J 

Code error estimate (sup norm 


. / -1.8 X 10-'' -6.9 X 10-« \ 
* -6.9 X 10-^ -2.6 X 10-^ ^ 
): 2.4 X 10-* 




Determinant of residue matrix: —7.2 x 10 — i2.0 x 10 




Value of |ai|: 1.3 x IQ-^. 
Integration tolerance for ^f: 10 


-9 




T. X or r^n^n.N / -297.2883 -110.8968 ^ 
Res(M,, A = 35.560604) = J 

Code error estimate (sup norm 


. / -5.7 X 10-^ -2.1 X 10-^ ^ 
+ * 1^ -2.1 X 10-^ -7.9 X 10-^ ^ 
): 4.0 X lO-"* 




Determinant of residue matrix: —7.3 x 10 — il.8 x 10 




Value of |ai|: 1.3 x 10"^. 
Integration tolerance for ^f: 10 


-9 





Table 3: Results for Equation 3 



are not known in closed form. However it is a relatively straightforward matter to compute 
approximations using the code SLEUTH |Q. For example, computing at different tolerances 
and using different truncations of [0, oo), the following approximations seem to be correct to 
all decimal places quoted: 

Ao = 35.560604, Ai = 128.113477, A2 = 297.84692. 

For Neumann boundary conditions y"'{0) = = y"{0) the corresponding approximate eigen- 
values obtained were 

Ao = 6.199245, Ai = 43.002631, A2 = 136.295990. 

^From this numerical evidence there is no overlap between the first few eigenvalues of the 
Dirichlet and Neumann spectra. To investigate whether or not there are likely to be HELP 
inequalities associated with this equation, we must compute the residues of the Titchmarsh- 
Weyl matrices at these eigenvalues using our code. The results are shown in Table |^. These 
residue matrices appear (to within the error we expected at the given tolerance) to be of rank 
1. This suggests that there is no HELP inequality associated with (|7^). 
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